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^ • ABSTRACT 

! A method based on Lucy (1974) iterative algorithm is developed to invert the equation 

of stellar statistics for the Galactic bulge and is then applied to the if-band star counts 
from the Two-Micron Galactic Survey in a number of off-plane regions (10° > |6| > 2°, 
\l\ < 15°). 

The top end of the K-hand luminosity function is derived and the morphology 
^ ' of the stellar density function is fitted to triaxial ellipsoids, assuming a non- variable 

luminosity function within the bulge. The results, which have already been outlined 
by Lopez-Corredoira et al. (1997b), are shown in this paper with a full explanation 
of the steps of the inversion: the luminosity function shows a sharp decrease brighter 
than Mk ~ —8.0 mag when compared with the disc population; the bulge fits triaxial 
' ellipsoids with the major axis in the Galactic plane at an angle with the line of sight 

Q>^ , to the Galactic centre of 12° in the first quadrant; the axial ratios are 1:0.54:0.33, and 

the distance of the Sun from the centre of the triaxial ellipsoid is 7860 pc. 

The major-minor axial ratio of the ellipsoids is found not to be constant, the best 
fit to the gradient being = (8.4 ± 1.7) x exp (-t/(2000 ± 920) pc), where t is the 
O ' distance along the major axis of the ellipsoid in parsecs. However, the interpretation 

| H I of this is controversial. An eccentricity of the true density-ellipsoid gradient and a 

C/3 . population gradient are two possible explanations. 

' The best fit for the stellar density, for 1300 pc < t < 3000 pc, are calculated for 

both cases, assuming an ellipsoidal distribution with constant axial ratios, and when 
■ Kz is allowed to vary. From these, the total number of bulge stars is ^ 3 x 10^" or 

' 4 X 10^°, respectively. 

b : 

5t 1 Key words: Galaxy: structure — infrared: stars — Galaxy: stellar content — stellar 

statistics 



1 INTRODUCTION 

This paper examines two aspects of the bulge: the luminosity 
function for the brightest stars in the K (2.2 ^m) band and 
the density distribution of these stars. 

There are many aspects of the bulge of the Galaxy that 
are still unknown, mainly because of the high extinction 
due to interstellar gas and dust. One of these unknowns is 
the near-infrared luminosity function, which has been princi- 
pally derived from observations in Baade's Window (Frogel 
& Whitford 1987; Davidge 1991; De Poy et al. 1993; Ruelas- 
Mayorga & Noriega- Mendoza 1995; Tiede et al. 1995). Gould 
(1997) and Holtzman et al. (1998) have used the Hubble 
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Space Telescope to study the V and I luminosity functions. 
However, extrapolations from Baade's and other clear win- 
dows to the whole bulge may not be appropriate, in partic- 
ular because these are 'special' regions. Furthermore, these 
regions are very small, containing relatively few stars, so 
they give very poor statistics at the brighter magnitudes. 
This bright end of the luminosity function is very important 
in order to determine the age of the population, for instance. 

Many authors have found non-axisymmetry in the 
Galactic bulge^ (Feast & Whitelock 1990) through the anal- 
ysis of star counts (Nakada et al. 1991; Weinberg 1992; 

t Some authors call it the "bar" instead of the bulge. See, for 
instance, Gerhard, Binney & Zhao (1998). 
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Whitelock et al. 1991; Stanek et al. 1994, 1996; Wozniak 
& Stanek 1996) or integrated flux maps (Blitz & Spergel 
1991; Weiland et al. 1994; Dwek et al. 1995; Sevenster 1996). 
This asymmetry has a negligible out-of-plane tilt (Weiland 
et al. 1994) and gives more counts in positive than in neg- 
ative galactic longitudes. However, other authors (Ibata & 
Gilmore 1995; Minniti 1996) claim that axisymmetry is suit- 
able. Besides the discussion about whether there is triaxial- 
ity or not, the actual shape and inclination of the bulge is 
also under debate with currently no clear agreement among 
different authors. 

Traditionally, star counts have been interpreted by fit- 
ting parameters to the functions. An assumption of an a 
priori shape of the bulge is made, along with the character- 
istics of its population. Free parameters are then fitted to the 
data and the model is obtained. This is the usual way of ex- 
tracting information concerning the different components of 
the Galaxy from star counts (Bahcall & Soneira 1980; Buser 
& Kaeser 1983; Prichet 1983; Gilmore 1984; Robin & Creze 
1986; Ruelas-Mayorga 1991; Wainscoat et al. 1992; Ortiz & 
Lepine 1993). The number of possible parameters to fit is 
limited to a priori assumptions about the shape (ellipsoidal, 
etc.) needed. 

In general, surface brightness maps are also interpreted 
by fitting parameters (Dwek et al. 1995; Freudenreich 1998). 
However, although these maps cover large brief ex- 

amination of the equations shows that they give no infor- 
mation on the luminosity functions. Therefore, when mak- 
ing the fit to the bulge, the number of free parameters is 
very small and applies only to the density function. Even in 
Freudenreich (1998), where in total some 30 parameters are 
used, only a very few of these apply to the bulge and only a 
very few parameters are solved at one go. 

In this paper we examine the TMGS star counts be- 
tween mK=4 and 9 mag in 71 regions across the bulge. 
The counts for these regions are shown in Hammersley et 
al. (1999), where a qualitative discussion on the counts is 
presented. It is shown that the counts are highly asymmet- 
ric in longitude when compared with the predictions of a 
symmetric model. 

Clearly, there is a relation between surface brightness 
and star counts as one is the integral of the other, but they 
are not the same or even similar and cannot be handled in 
the same manner. One way of looking at the difference be- 
tween the two is to consider that at a single position a surface 
brightness map gives just a single value whilst star counts 
give a counts, vs magnitude plot. From this plot alone it is 
possible to determine things about the structure. Therefore, 
while star counts and surface brightness maps are clearly 
related, they behave very differently. Many authors, includ- 
ing ourselves, have already used the fitting approach to look 
at the surface brightness COBE-DIRBF, data (we note that 
Binney et al. 1997 have tried inversions on the surface bright- 
ness maps), but this is not the best for star counts in the 
present situation. 

One of the major advantages in analysing star counts as 
opposed to surface brightness maps is that the magnitude 
range can be limited in order to highlight the features of 
interest. This is of particular value when looking for triaxi- 
ality because if one region is significantly closer than another 
then, simply from the inverse square dependence with the 
distance, the sources from the further region are not de- 



tected until a fainter magnitude. Hammersley et al. (1999) 
show that in the TMGS star counts the size of the asym- 
metry amounts to some 50% of the bulge counts in some 
magnitude ranges, in the CO-B-&DIRBE maps the asymme- 
try is far less. For this reason analysis of 2-/xm star counts in 
a certain magnitude range will be far more sensitive in de- 
termining the triaxiality of the bulge than surface brightness 
maps. 

Whilst large-area star counts, as used here, contain far 
more information than surface brightness maps they are in- 
trinsically far more difficult to analyse. A priori, neither is 
known and furthermore there is no reason to believe that the 
luminosity function (LF) is a simple analytical expression. 
Therefore, whereas fitting a surface brightness map there 
will only be a few free parameters, this is not the case for 
star counts. In this case the number of free parameters would 
rise unmanageably and so we would be forced to adopt a pri- 
ori assumptions on the LF and density functions with the 
severe risk that the final result is dependent on these initial 
assumptions. 

We have therefore chosen a different approach, that of 
direct inversion. Assumptions on the shape of the solution- 
functions (in this case, these are the luminosity function and 
the density of stars) are not made but instead come directly 
from the data by means of an "inversion" technique. Once 
the solutions for the functions are produced by the inversion, 
they are compared a posteriori to some known analytical 
expression (for instance, an ellipsoidal shape for the bulge 
isodensity contours) and, afterwards, fitted to them. This 
method allows all possible solutions to be examined, rather 
than solely that of the initial assumption and the only fitting 
are density contours to a density map. No attempt need to 
be made to fit a density function to the star counts. 

Since the first decades of this century, attempts have 
been made to invert the star-count equation (eq. |^). How- 
ever, problems such as excessive patchiness of extinction in 
optical star counts or instabilities of an ill-posed problem 
in the mathematical technique of inversion, hindered the 
development of the technique. In this paper the extinction 
problems are ameliorated by using the near-infrared K band 
and the instabilities by using a statistical iterative algorithm 
of inversion (Lucy 1974) . A full explanation of the inversion 
is developed in this paper (with the core in ^ whose re- 
sults have already been outlined by Lopez-Corredoira et al. 
(1997b). 



2 NEAR-INFRARED DATA 

K-hand star counts were taken from the Two Micron Galac- 
tic Survey (TMGS; Garzon et al. 1993, 1996), which covers 
about 350 deg^ of sky and has detected some 700000 stars 
in or near the Galactic plane. This survey provides if-band 
observations of several regions that cross the Galactic plane, 
in the areas -5° < / < 35°, |&| < 15° and 35° < / < 180°, 
l&l < 5°. 

Regions from three strips of constant declination are 
used (Table |l]) In this study of the bulge. More specif- 
ically, 71 regions were selected from those strips in off- 
plane regions, but not too far from the Galactic centre 
(10° > |6| > 2°, \l\ < 15°). Each region has an area on the 
sky between 0.4 and 1.9 deg^. The chosen regions are listed 
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Table 1. Constant-declination TMGS strip used in this paper. 

'5ccntral(J2000) Cut in the Galactic plane Strip width (A<5) 
(deg) (deg) 

-29°43'32" / = -0.9 2.51 

-22°26'40" 1 = 7.5 1.63 

-15°33'24" I = 15.4 0.78 



l=-2.5 deg., b=2.1 deg. 




Figure 1. Comparison of cumulated star counts without con- 
fusion correction and those corrected according to the method 
explained in Lopez-Corredoira et al. (1997a) with a linear ex- 
trapolation of the differential star counts over magnitude 9.4. 



in Table ^. There is an overlap in the neighbouring regions 
such that some stars fall into two regions. The total covered 
area of sky covered is 75 deg^. This area is far greater than 
that used in Baade's window or any of the other low extinc- 
tion region and hence provides much better statistics for the 
top end of the bulge LF. 

The chosen regions contain principally bulge and disc 
stars. The area near the Galactic plane was not used in 
order to avoid components which belong neither to the bulge 
nor to the disc (e.g. spiral arms) and the high and variable 
extinction. The outer limits were set so that the bulge-to- 
disc stellar ratio was still acceptable, i.e. so that there were 
sufficient bulge stars in comparison with disc stars to make 
the study of the bugle meaningful.. 

The survey is complete between the magnitude limits 
rriK = 4.0 mag and rriK ~ 9.2 mag, except for the regions 
very near the Galactic centre where source confusion re- 
duced the faint limit by about half a magnitude, although 
the detection limiting magnitude of the survey is in excess 
of 10 mag. Hence, inversion will be applied up to rriK = 8.6 
mag for the regions of the strip with declination —30° and 
up to rriK = 9.0 mag for the remaining cases. 

Figure ^ shows cumulative star counts, A^, for the three 
strips up to rriK = 9 mag as a function of b {I also varies, as 
can be seen in Fig. 

Within this range of magnitudes, confusion effects are 
negligible. This was determined from the application of the 
method explained by Lopez-Corredoira et al. (1997a) by as- 
suming an extrapolation to fainter magnitudes (Fig. hh. As 



can be seen in the figure, the counts are nearly the same with 
or without correction. That confusion is not significant for 
the areas chosen can also be seen in the figures in Hammer- 
sley et al. (1999) where the TMGS star counts are directly 
compared with the W92 model counts (Cohen 1994). Taking 
into account that the correction is based on an extrapola- 
tion and the changes are minor when compared to the other 
sources of error, it is preferable to avoid any correction and 
use the original counts. 



3 THE STELLAR STATISTICS EQUATION 
FOR THE BULGE 

3.1 Cumulative star counts 

For each of the 71 regions centred on galactic coordinates 
{I, b)i, where i is the field number, the cumulative star counts 
observed in a filter K, Nk, up to a magnitude rriK in a given 
region of solid angle uj is the sum of the stars over the beam 
with such an apparent magnitude (Bahcall 1986). Assuming 
a luminosity function which does not vary with the spatial 
position for each Galactic component c, this is 



/•oo 

„ Jo 



a Kir)) 



X Dc{r)r^dr, 
where 



Mk 



^kAMk)= / <t)KAM)dM.., 



(1) 



(2) 



(j}K,c is the normalized luminosity function for the K band 
in the component c; Dc is the density function in the com- 
ponent c and aK{r) is the extinction along the line of sight 
for the K band. 



3.2 Extinction 

If the star counts, Nk, for eq. (^), and the luminosity func- 
tion, 4iK,c, are known then the densities and the extinction 
would be the unknown functions. The extinction can be sep- 
arated from the last integral equation by means of a suitable 
change of variable (Bok 1937; Trumpler & Weaver 1953; Mi- 
halas & Binney 1981, ch. 4): 



t r\0.2a (r ) 

PK = 10 'r. 



:[PA-(r)] = 



D.{r) 



[l + 0.2(ln lQ)r!i2^P-) IQO.aaK^r) ' 
which transforms the stellar statistics equation into 



(3) 
(4) 
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Table 2. The regions whose star counts are used to invert and extract information about the bulge. 



I 

(deg) 






; 

(deg) 


(deg) 




(deg) 


(deg) 




-6.3 


7.8 


0.4 


2.6 


-6.1 


1.9 


9.0 


-2.7 


1.4 


-5.7 


7.1 


0.8 


3.0 


-6.9 


1.9 


9.1 


-2.9 


1.4 


-5.2 


6.4 


1.3 


3.4 


-7.7 


1.8 


9.2 


-3.0 


1.4 


-4.7 


5.7 


1.8 


3.8 


-8.5 


1.8 


9.6 


-3.9 


1.4 


-4.2 


5.0 


1.9 


4.2 


-9.2 


1.8 


10.1 


-4.7 


1.4 


-3.7 


4.3 


1.9 


1.3 


9.9 


1.4 


10.5 


-5.5 


1.4 


-3.2 


3.5 


1.9 


1.8 


9.1 


1.4 


10.9 


-6.3 


1.4 


-2.7 


2.8 


1.9 


2.3 


8.4 


1.4 


11.3 


-7.1 


1.4 


-2.6 


2.7 


1.9 


2.9 


7.6 


1.4 


11.7 


-8.0 


1.4 


-2.5 


2.5 


1.9 


3.4 


6.8 


1.4 


12.2 


-8.8 


1.4 


-2.4 


2.4 


1.9 


3.9 


6.0 


1.4 


12.6 


-9.6 


1.4 


-2.3 


2.2 


1.9 


4.4 


5.3 


1.4 


9.7 


9.9 


0.7 


-2.3 


2.1 


1.9 


4.9 


4.5 


1.4 


10.2 


9.1 


0.7 


0.3 


-2.0 


1.9 


5.4 


3.7 


1.4 


10.7 


8.3 


0.7 


0.4 


-2.2 


1.9 


5.8 


2.9 


1.4 


11.2 


7.4 


0.7 


0.5 


-2.3 


1.9 


5.9 


2.7 


1.4 


11.7 


6.6 


0.7 


0.6 


-2.5 


1.9 


6.0 


2.6 


1.4 


12.2 


5.8 


0.7 


0.7 


-2.6 


1.9 


6.1 


2.4 


1.4 


12.6 


4.9 


0.7 


0.7 


-2.8 


1.9 


6.2 


2.3 


1.4 


13.1 


4.1 


0.7 


0.8 


-2.9 


1.9 


6.3 


2.1 


1.4 


13.6 


3.3 


0.7 


0.9 


-3.1 


1.9 


8.7 


-2.1 


1.4 


14.1 


2.4 


0.7 


1.3 


-3.9 


1.9 


8.8 


-2.2 


1.4 


14.2 


2.2 


0.7 


1.8 


-4.6 


1.9 


8.9 


-2.4 


1.4 


14.3 


2.1 


0.7 


2.2 


-5.4 


1.9 


8.9 


-2.6 


1.4 









X Ac,K{pK)pKdpK- (5) 

The functions Ac,k{pk) do not have a direct physical 
meaning but are fictitious densities as a function of a fic- 
titious distance which coincides with the real distance only 
when there is no extinction (see Calbet et al. 1995). 

For the extinction we have followed Wainscoat et al. 
(1992, hereafter W92), who assume that the extinction has 
an exponential distribution with the same scale length as 
the old disc, 3.5 kpc, and a scale height of 100 pc. This is 
normalized to give Ak = daj-c/dr = 0.07 mag kpc~^ in the 
solar neighbourhood. Although this model is crude it is suf- 
ficient for our purposes. As the areas of interest are off the 
plane, the extinction in the direction of the bulge sources 
is between 0.05 to 0.5 mag at K (ten times lower than in 
V). The evidence from the 2.2-/im surface brightness maps 
is that there are off-plane clouds, but these are isolated so 
if a strip did cross a cloud it would affect only one or two 
regions which would have a minor effect on the final result. 
In fact, Hammersley et al. (1999) show that in the regions 
chosen there are no major dips in the counts and hence no 
isolated clouds. Furthermore, in this paper there is a dis- 
cussion on the IR extinction in the plane and comparison is 
made with the W92 model, which uses the above model for 
the extinction. It is shown that in the solar neighbourhood 
this model works well and remains valid to a galactocentric 
distance of about 4 kpc where the molecular ring is situated. 
Inside the ring the extinction is then over estimated. How- 
ever, it should be noted that for the lines of sight used here 
the majority of the extinction occurs in the first few kpc, i.e. 
while the line of sight is close to the Galactic plane. There- 



fore, the extra extinction added by the model in the inner 
galaxy is a small proportion of the total extinction along 
the line of sight, which is in turn already small. This effect 
is clearly demonstrated by Hammersley et al (1999) for the 
I — 7° strip where the effect of the overestimated extinction 
can only be seen within 0.5° of the plane. 

Another possible cause for concern could be if there was 
a general asymmetry in the extinction, either from above to 
below the plane or between positive and negative longitudes. 
However, it must be noted that the analysis of Freuden- 
rich (1998) of the COB£^DIRBE surface brightness maps 
shows no such asymmetry. Furthermore, Hammersley et al. 
(1999) have analysed the asymmetry in the TMGS bulge 
star counts and show that the form is not consistent with the 
asymmetry being caused by extinction. Therefore, although 
the extinction model is crude it is valid for the purpose used 
here. 

So, from aK, the relationship is obtained between A 
— the fictitious density — and D — the real density — for 
each component, using eq. (|^; therefore eq. (^ will be used 
hereafter. 



3.3 Subtraction of the disc 

The components cannot all be solved simultaneously and 
the inversion of eq. (^ can only be solved when the number 
of components, c, is restricted to one. 

It will be assumed that, in the chosen regions, the con- 
tribution to the star counts will be primarily from the disc 
and bulge. In order to isolate the bulge component, there- 
fore, the contribution of the disc must be subtracted from 
the total counts for each region. 

The model of the disc coded by us was based on W92, 
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Figure 2. N{mK = 9.0 mag) along the three strips that are used with constant declinations: 5 = —30°, which cuts the plane at Z = — 1°; 
5 = —22°, which cuts the plane at Z = 7°; and S = —16°, which cuts the plane at / = 15°. 



which follows Bahcall & Soneira (1980). It has been used 
because it provides a good fit to the TMGS counts in the 
region where the disc dominates (Cohen 1994b; Hammersley 
et al. 1999). The W92 model was revised by Cohen (1994a) 
but this does not significantly alter the form of the disc in 
the areas of interest. Three examples of those fits are shown 
in Fig. in regions where the disc is isolated (note that 
the regions used in these plots are different from the regions 
used for the inversion specified in §^). A more detailed com- 
parison of the W92 model and the TMGS is presented by 
Hammersley et al. (1999), who examine some 300 square 
degrees of sky. Hence, by extrapolation, it is expected that 
this disc model will adequately reflect the disc components 
along the lines of sight used in this paper. Initially, it was 
also expected that the W92 model would give an adequate 
fit for the bulge counts; this, however, was not the case as 
can be clearly seen in Hammersley et al. (1999). 



3.4 Fredholm integral equations of the first kind 

Once the disc star counts are subtracted, a Fredholm integral 
equation of the first kind is derived (see Trumpler & Weaver 
1953): 

NK,buigc{mK) = NiiimK) - NK,disc{mK) 



= a; / «&x,buigG(»nK + 5 — Slog^o PJf)^buigc,if (Pi"f) 

X PKdpK, (6) 

where A is the unknown function and $ is the kernel of the 
integral equation. 

When "1> is the unknown function instead of A, then 
a new change of variable can be made: Mk ~ rriK + 5 — 
Slogj^Q pK, and a new Fredholm equation of the first kind is 
obtained: 

NKimx) = 200(ln 10)10^ 

Ak(10 B )10—^ $K{MK)dMK. (7) 

- oo 

In this case, the kernel is A instead of $. Any method of 
inverting eq. (^) is also applicable to this integral equation 

(§)• 

4 INVERSION OF THE STELLAR STATISTICS 
EQUATION 

The inversion of integral equations such as (^) or ^ is ill- 
conditioned. Typical analytical methods for solving these 
equations (see Balazs 1995) cannot achieve a good solution 
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iV^ ^5^.3630 

Area'^ 1.00 sq deg 
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K Ma£>nitude 

a) 




12 11 10 9 8 7 6 5 4 

K Maenitude 

c) 



Figure 3. Differential star counts, the derivative of the cumu- 
lative star counts. Rhombi are TMGS data. Lines represent the 
W92 model: the solid line stands for counts for all components; 
the dotted line stands for disc counts; long-dashed line for spiral 
arms; short-dashed and dotted line for the ring; shot-dashed line 
for the bulge; long-dashed and dotted line for the halo. In these 
cases - a), b) and c) - disc and total counts are nearly coincident 
because the disc gives the most part of the stars. 




b = -B-OCr 
1 J 24.119 
Area = 1.00 "^qNleg 
Solar displacement^^x 15.00 po.^ .-■ 



7 

Magnitude 



b) 



because of the sensitivity of the kernel to the the noise of 
the counts (see, for instance, Craig & Brown 1986, ch. 5). 

Since the functions in these equations have a stochastic 
rather than analytical interpretation, it is to be expected 
that statistical inversion algorithms will be more robust. 
This is confirmed by several authors, for instance Turchin 
et al. (1971), Jupp et al. (1975), Balazs (1995). 

From among these statistical methods, we have selected 
Lucy's algorithm (Lucy 1974; Turchin et al. 1971; Balazs 
1995), an iterative method, the key to which is the inter- 
pretation of the kernel as a conditioned probability and the 
application of Bayes' theorem^. 

t Bayesian methods have multiple applications in astrophysics. 
Inversion problems are particular cases of these applications 
(Loredo 1990). 



In eq. (^ , A is the unknown function, and the kernel is 
$, which depends on the apparent magnitude conditioned to 
the fictitious distance p. The fictitious density A can also be 
understood in terms of a probability density (the probabilit; 
of finding a star with fictitious distance p). Thus, eq 
can be rewritten as (hereafter, the notation for component 
or passband will be dropped) 



N(m) 



where 



Aip)P{m\p)dp, 



P{m\p) = p $(m + 5-5 logn, p)- 



(8) 



(9) 



The inverse conditioned probability, i.e. the probability of 
star being at a fictitious distance p, once its apparent mag- 
nitude m is known, is given by Bayes' theorem: 



Q{p\m) = 



A(p)P(mlp) 



A{x)P{m\x)dx' 
From the definition of conditioned probability, 
A(p)P(m!p) = N(m)Q{p\m), 
and, hence, we get directly: 

J™-^" dmN{m)Q{p\m) 



Hp) 



XT""" dmP{m\p) 



(10) 



(11) 



(12) 



Equations ( p^ ) and ([To|) together lead to an iterative 
method^ of obtaining the unknown function A(p): 



'(p) = A'-(p) 



rn 



N° 



^P(m|p)dm 



(13) 



where A'^"'''^ represents the observed cumulative counts and 



iV(m) 



/\''(x)P(m\x)dx. 



(14) 



5 For the numerical calculation of these integrals p is placed into 
discrete logarithmic intervals (the (m, log tt) method; Mihalas & 
Binney 1981, ch. 4) in such a way that log^o Pjf is regularly 
spaced. 
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This development is more general than Lucy's. Lucy's 
algorithm (Lucy 1974) 's algorithm was expressed for cases 
with /J^""" P{i^\p) = 1, which is not true in the case dis- 
cussed here because the range of magnitudes is limited. The 
need for the denominator in eq. ( [l^ ) was already recognized 
by Scoville et al. (1983). 

The iteration converges when iV = iV°''=, i.e. when 
A''^^ = A"". The first iterations produce a result which is 
close to the final answer, with the subsequent iterations giv- 
ing only small corrections. 

This algorithm has a number of good properties (Lucy 
1974, 1994): both the luminosity function and the density 
are defined as being positive, the likelihood increases with 
the number of iterations, the method is insensitive to high 
frequency noise in A'^"'''', etc. 
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4.1 stopping criteria for the iterative process and 
initial trial solution 

From Lucy (1994), the appropriate moment at which to stop 
this kind of iterative process is when the curvature of the 
trajectory in the H-S diagram is a minimum. H and entropy 
(S) are defined by: 



i 



and 

s = -^Ann§, 



respectively, and the curvature in the H-S diagram is: 
\S'H" -H'S"\ 



(5-/2 + ///2)3/2 



(15) 



(16) 



(17) 



where the derivatives are with respect to the number of it- 
erations, and the sums over i and j correspond to discrete 
values of the p and m integrals respectively. 

Tests were carried out on the data set using this crite- 
rion (see an example in Fig. In general there is a mini- 
mum after three iterations, corresponding to a non-relaxed 
state of the process. Afterwards, k, is increase up to around 
10 iterations, where it then falls off again to a minimum, 
and then increases again. Apart from first minimum at 3 
iterations, the most relevant minimum seems to be that at 
around 10 iterations. 

However, this criterion is not very accurate for the nois- 
iest cases and on occasions the last iteration may not be the 
most appropriate one to end at. Occasionally, it stops too 
early and therefore hinders the extraction of further infor- 
mation that could be exploited. 

Therefore, the following criteria are adopted for ending 
the iterations: 

(i) The number of iterations must be greater than 10 and 
smaller than 10000. The process will always be stopped when 
the number of iterations exceeds 10000. The A' variations 
are too small after 10000 iterations, so no more are made. 

(ii) For fewer than 1000 iterations, the iterative process is 
stopped when the solution is within the noise, i.e. when the 
average over m of the distance between N^{m) and N°^^{m) 



Figure 4. Curvature versus iteration number in an inversion case. 



is less than the average over m of a random noise with Gaus- 
sian distribution of N°'°^{Tn) with am = S{N°'°^{m)), the 
Poissonian noise of N°^^{m). 

This last point will be clarified and the numerical algo- 
rithm to be used explained in what follows. N[ (the subindex 
i stands for the discrete value of m) is at Si as from N°^^, 
i.e. 



Si 



(18) 



The normalized probability of a point at distance Si as 
from its real value is 



Piisi) = erf(si). 



(19) 



where erf(2;) — {2/^/tt) e~'" du is the error function. 
Thus, since the Pi distribution is nearly uniform between 
and 1, then the Si distribution follows 



Pidpi = -. 



(20) 



'Nearly' because it is exact when n oo, and there are 
some fiuctuations when n is not too large. 
Thus, within the noise means that 



n 6 



(21) 



and this is second stopping criterion. 

The sum of pi is calculated instead of the sum of Pi be- 
cause the difi'erence distribution is not exactly Gaussian and 
a power of pi gives a higher weighting to the large deviations 
(larger than 1-2 a) . In any case, this is only an approximate 
criterion. 

The final solution does not depend on the initial trial 
solution, A*'^, when the number iterations is high enough. 
However, A^"^ may approach A'^"'''' in a different way depend- 
ing on the initial trial solution when the noise of the counts 
is high, because the process is stopped after a few iterations 
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which will give slightly different solutions. In order to avoid 
this inffuence for the noisiest data to be invertedp], the trial 
solution was fed back with the smoothed result of the pre- 
vious inversion and inverted again. As will be discussed in 



4, three inversions are made. In the second and the third 
iterations the trial solutions are fed back with the previ- 
ous outcome, once it has been fitted to a smooth analytical 
function (in this case ellipsoids, as seen in ^3). 



4.2 Distance range 

As the case described here is the application of the method 
to the bulge of our Galaxy, the numerical calculation of the 
distance integral are carried out over 2000 pc< pK < 30000 
pc, as all of the stars are known to be contained within this 
distance. The real distance, r, is somewhat lower than p (see 
eq. P), but the difference is small for low- extinct ion regions 
such as those used here. 

It noted that, following numerical experiments with 
Lucy's algorithm, the minimum distance has to be kept 
within tolerable limits. Spurious fluctuations arise when 
small distances are included. This is related to the propor- 
tionality of the kernel to p^, so that large variations in the 
density at small distances do not significantly change the 
number of counts. 

A similar problem arise for the maximum distance to 
which the sources can be distributed. If the maximum limit 
is too large then a spuriously high density might appear at 
large distances. The reason is that very distant stars should 
be very luminous to be observed and, since the luminosity 
function for very luminous stars is very small, any sources 
placed at a large distance will lead to a high density at that 
distance. 

The application of this method to the bulge does not 
lead to problems since the distance range is known to be 
limited: the Galaxy has a boundary and the number of bulge 
stars in the solar neighbourhood is negligible. Nevertheless, 
it should be noted that care should be taken before applying 
this method to other Galactic components. For instance, it 
is possible that inverting the counts to obtain the Galactic 
disc density could encounter the above problems. 
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Figure 5. Recovery of the theoretical luminosity function 
through the inversion process. Three cases: a), b), c). 
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4.3 Example of application 

Inversion of the stellar statistics equation has been discussed 
by many authors, much more often in theory than in prac- 
tice, and doubt has been cast on the viability of such an 
inversion. It has even even said that as the solution is non- 
unique (Gilmore 1989), which would lead to instability in 
the inversion. Except for some particular kernel functions 
(Craig & Brown 1986; ch. 4) this is not in fact the case, 
as we shall attempt to demonstrate here. The question of 
uniqueness is important only from a theoretical standpoint. 
In practice, the only relevant issue is whether the method 
is able to obtain a solution close to the real one when the 
counts are affected by noise, which always produces devi- 
ations from the real solution. The important thing is that 

^ Very noisy data are eliminated. In this case, the 37 least 
noisy regions out of 71 are used when the density is the unknown 
function. 



this solution be not very far from the true solution. That 
the solution is not unique need not be important when all 
solutions be close to each other. 

In order to test the reliability of the method, a num- 
ber of simulations were made. A luminosity function and a 
fictitious density function were constructed. The cumulative 
count per square degree, N{m), were then calculated by in- 
tegrating eq. A random noise with a Gaussian distribu- 
tion is added to each bin. The cumulative counts with noise 
are then represented by N°^^{m). When Lucy's algorithm is 
applied, with the same luminosity function and A^(p) = 1 
(the choice of the trial initial solution does not affect the 
outcome), the results shown in Fig. ^ are obtained. 

The inversion is not perfect since it is affected by the 
noise, but the results are fairly good. There is a "hump" 
in the first case at short distances, and a large increase in 
the density at large distances in the second case. However, 
it should be noted that the hump is at the 10% level of 



© 0000 RAS, MNRAS 000, 000-000 



Inversion ... Bulge 9 




c) 



the primary peak, which is located very close to the correct 
distance of 10000 pc. Similarly, the excess at large distances 
is at the level of only a few percent of the peak sources. 
Sensitivity to noise is higher for distances less than 7000 pc 
or greater than 15000 pc, as explained in §^^, and this is 
reproduced in the experiments. 

When the method is asked to recover two peaks, the 
inversion gives poorer results. However, were the number of 
iterations increased beyond the 10000 limit, then the second 
peak in Fig. |^ c) would rise and become closer to the orig- 
inal. Again, however, both peaks are correctly located and 
the total number of sources in each peaks is very close to 
the original. Apart from these details, the general shape of 
the peaks is recovered. Other numerical experiments were 
performed with similar results. 

The bulge is a single-peaked structure so the proposed 
stopping criteria are sufficient. Since noise is random, the 
composition of the three-dimensional densities from the in- 
version for different regions (/,&) will attenuate the average 
deviations. 

Application to equation (^) , instead of ^ , deserves sim- 
ilar considerations. 



4.4 Method of deriving both the luminosity 
function and the density 

The equations and (Q) can be solved for either the lu- 
minosity function or the density function, but not for both 
simultaneously for each region. Since both functions, A and 
are of interest but accurate information is not available 
for either of them, the following method was used. 

To begin with, a first order approximation for the den- 
sity was assumed. It was taken from the axisymmetric W92 
model. A simple comparison showed that the W92 luminos- 
ity functions suggested that there were possible problems 
with the brightest sources, although the density function 
did give a reasonable starting point. Therefore, it was de- 
cided to solve first for the average luminosity function using 
the W92 bulge density. 

With this density distribution, eq. (Q) is inverted by 
means of Lucy's algorithm to provide the luminosity func- 



tion for each of the regions {I, b) in Table g. The weighted 
average of all luminosity functions was then calculated. 

We have made the assumption that the bulge luminosity 
function is independent of position. This assumption is sus- 
pect (see Frogel 1988, Section 3) since the observed metal- 
licity gradient might affect the luminosity of the AGB stars, 
although not the non-variable M-giants whose bolometric 
luminosity function is nearly independent of the latitude 
(Frogel et al. 1990). Some authors claim that there is a pop- 
ulation gradient (Frogel 1990; Houdashelt 1996; Frogel et al. 
1999), while others do not (Tyson & Rich, 1993, show that 
there is no metallicity gradient up to 10° out of the plane; 
Ibata & Gilmore, 1995, argue that there is no detectable 
abundance gradient in the Galactic bulge over the galacto- 
centric range from 500 to 3500 pc). While the assumption 
may not be strictly true, it is nevertheless a useful approxi- 
mation in deriving mean properties of the bulge. 

With this averaged luminosity function, eq. (P) was in- 
verted to derive a new density distribution by means of 
Lucy's algorithm for each region. In this step the 37 regions 
with the highest counts were used, as the determination of 
the density is more sensitive to noise. 

The inversion of the luminosity function is more stable 
because the density distribution is sharply peaked and so 
the kernel in eq. (M) behaves almost eis a Dirac delta func- 
tion. Hence, the shape of the density distribution does not 
significantly affect the shape of the luminosity function. 

The new density was then used to improve the lumi- 
nosity function, etc. The whole process was iterated three 
times, which was enough for the results to stabilize as can be 
seen in Fig. ^: it is seen how the result of the third iteration 
is very close to the first, i.e. stabilization is reached in the 
first iterations. This small variation in successive iterations 
i s real ly a convergence to the solution since, as is shown in 
^.3.4 and 8:6.4.1, the counts are approximately recovered 
when we project the bulge obtained from the inversion. 

The functions of interest are (j), the derivative of $, and 
D, related to A by the change of variable expressed in eqs. 
(I) and (I). 



5 THE TOP END OF THE K LUMINOSITY 
FUNCTION 

After three iterations the luminosity function was nearly 
independent of the position {l,b)i, stable and hardly changed 
from the solution of the second iteration. Compare the first 
three iterations in Figure |^. In fact, even the first iteration 
came close to the final solution. 

The obtained luminosity function is shown in Fig. |^ and 
in table|| The derivative, 4>, of ${Mk), from eqs. (^ and (^ 
is the normalized probability of having absolute magnitude 
Mk per unit absolute magnitude. 

Figure ^ shows that for —10 mag < Mk < —8 mag the 
bulge luminosity function is significantly lower than that of 
the disc (Eaton et al. 1984). Hence, the density of very bright 
stars in the bulge is much less than in the disc. Fainter than 
Mk ~ —8 mag the luminosity functions of the disc and the 
bulge coincide, in agreement with Gould (1997). The lumi- 
nosity function for -10 mag < Mk < —8 mag is significantly 
below the synthesized luminosity function assumed by W92 
for the bulge in their model of the Galaxy (this can also be 
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Table 3. if-band luminosity function for bulge stars. 
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Figure 6. Luminosity function in the first three iterations. 



clearly seen in the W92 model and the TMGS in Hammer- 
sley et al (1999)). This discrepancy could arise from their 
not having taken into account that the brightest stars in the 
bulge are up to 2 mag fainter than the disc giants (Frogel & 
Whitford 1987). This would shift the W92 luminosity func- 
tion to the right in Figure ^ It should be remembered that 
the W92 model was developed to predict the IRAS source 
counts. IRAS could see only the very top end of the bulge 
luminosity function, and the sources responsible are all dust- 
shrouded AGB stars. The dust enormously brightens the 12 
and 25 micron fluxes over the expected photospheric flux. 
In fact, at the distance of the bulge, IRAS could not see 
purely photospheric stars at all. The TMGS, however, can 
detect normal bulge M giants (Frogel & Whitford 1987), 
not only AGBs, and the presence of dust leads only to a mi- 
nor increase in the K brightness. Therefore, in the TMGS 
while it is true that we do see the extreme AGB stars de- 
tected by IRAS, they in fact represent only a tiny fraction 
of the detected sources in each magnitude bin. Hence, the 
top end of the IRAS luminosity function and the top end of 
the TMGS luminosity function are dominated by different 
types of sources and so W92 could be close for IRAS but 
not get the top end of the K star counts correct. 

Between Mk = —8 mag and Mk = —6 mag (corre- 
sponding to the fainter limit of the TMGS at the distance of 
the bulge) the luminosity function of W92 does coincide with 
that determined here. As has already been noted, the result 
from the first iteration of the luminosity function (when the 
assumed density function was that of W92) is very close to 
the final result, particularly for absolute magnitudes fainter 
than Mk < —8 mag. This implies that for the lines of sight 
used here the W92 model does correctly predict the num- 
ber of bulge stars per magnitude per square degree for —8 
mag < Mk < —6 mag, even though this model was aimed 
at matching the IRAS source counts. Given the match over 
this magnitude range we have chosen to use the W92 lumi- 
nosity function for the magnitudes fainter than Mk ~ —G 
mag, so that the luminosity function can be normalized. 

Comparison with the bolometric luminosity function 
obtained by other authors (see references in the introduc- 



tion) is not possible since bolometric corrections are not 
available. Also, in most of cases the magnitude interval is 
different. Tiede et al. (1995) provide, by combining data 
from different works, the luminosity function in the K band 
as a function of the apparent magnitude in the range 5.5 mag 
< rriK < 16.5 mag. The brightest magnitudes are taken from 
Frogel & Whitford (1987). The comparison with our lumi- 
nosity function is not direct since they have not normalized 
their luminosity function to unity; moreover, they have not 
taken into account the narrow but non- negligible dispersion 
of distances. In Fig. 16 of Tiede et al. (1995) there is a fall- 
off in the luminosity function for rriK < 6.5 mag or in Fig. 
18 of Frogel & Whitford (1987) for Afboi < -4.2 mag, which 
could be comparable with that of our luminosity function at 
Mk ~ —8.0 mag. However, because of the much larger area 
covered by the TMGS, the error for the brightest magni- 
tudes is far lower in this paper, the result being pushed well 
above the noise; this is not the case for Frogel & Whitford 
(1987). 

The presented luminosity function for very bright stars 
(brighter than Mk ~ —9.5 mag) is of low precision. The 
number of bulge stars in this range is very small, so even 
small errors due to contamination from the spiral arms will 
mean that the luminosity function is overestimated and so 
the values should be taken as an upper limit. 



5.1 Age of the bulge 

The age of the bulge is an open topic. There are authors 
who think the bulge is older than the halo (Lee 1992) whilst 
others hold the opposite opinion (Rich 1993) . Although from 
the work presented here an accurate value for its age can- 
not be determined, the bulge is clearly older than the disc. 
The lack of very luminous stars in the bulge means that 
there are few supergiants and bright giants, and hence star 
formation regions. A comparison between the K-hand lumi- 
nosity function derived here and models of stellar evolution 
could provide some further clue in this controversial subject. 
The model of Bertelh et al. (1994), with a 10-Gyr popula- 
tion and solar metallicity, predicts that all the stars should 
be fainter than Mk = —8 mag, while these data show that 
there are some sources of —9.5 and —8 mag. This may indi- 
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Figure 7. Luminosity function in the K-hand (solid line). Com- Figure 8. Density along the line of sight in the region [l = 5.4°, 

parisons with W92 in the bulge and Eaton et al. (1984) in the b = 3.7°). 

disc are also provided. 



cate a mixture of populations with different ages embedded 
in the bulge. 



6 DENSITY DISTRIBUTION 

6.1 Density along the line of sight 

The second result is the density D{r) for each region {l,b), 
i.e. some points of the function D{r) — D{r,l,b). A is ob- 
tained by inversion of eq. (^) and then changing the variable 
in eq. (Bl) to recover D{r). 

As an example, the density distribution along the line 
of sight for one region (/ — 5.4°, b = 3.7°) is shown in 
Fig. H after extinction correction. As can be seen, the bulge 
distribution of stars has a maximum around 8 kpc. There 
is a rise from ~ 5 kpc to ~ 8 kpc, and a fall off after this. 
Similar results were obtained in the other regions, except 
for some fluctuations due to errors (the errors in the counts 
may provide this fluctuation; see §W). The 37 regions used 
were the least noisy and least affected by patchy extinction. 



6.2 Bulge cuts 

As was said in the regions come from strips with constant 
declination: 5 = -30°, 5 = -22° and S = -16°. The 37 
regions used for density inversion are come from the strips 
at S — —30°, S = —22° (as the bulge source density by 
S — —16° is low). A strip can be thought of as a surface in 
space (Fig. ^ one axis is in R.A. (i.e. constant declination) 
and the other is distance along the line of sight, which can be 
converted to a distance parallel to the Sun-Galactic center 
line. Figures ^ and ^ show these plots with the z-axis 
representing the density. Note that the density scale (height) 
is different in both figures. 

As can be seen, there are two peaks and a valley in both 
figures. The valley only indicates the absence of data due to 
the fact that the Galactic plane between b — —2° and b = 2° 
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Figure 9. Two constant-declination strips that cut the disc. The 
striped region is the Galactic plane zone, which was excluded. 



was avoided. If the plane data were included there would be 
only one peak. 

Galactic longitude increases and latitude decreases with 
increasing x. In Fig. the left side (negative x) of the 
valley has a lower density than the right side (positive x) 
due to the abrupt fall-off of the density with distance from 
the Galactic center. This is not observed in Fig. |l^ because 
this strip almost cuts across the Galactic center so both sides 
of the valley are nearly symmetric. 

When comparing the position of the peaks, and hence 
the maximum density, in both figures, the peaks are notice- 
ably closer to the Sun for S = -22° (l = 7.5°) than for 
5 — —30° {I — —1°). The non-axisymmetry of the bulge 
is the most plausible explanation for this and the bulge is 
closer to us at higher galactic longitudes. This can, in fact, 
be seen in the individual strips, as the left peak (i.e. larger 
I) is closer than the right one in both figures. Hammersley et 
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Figure 10. Plot of the density (height) as a function of both 
spatial coordinates defined by a cut of the bulge in 5 = —30°. 
Galactic latitude is increased from left to right (a;-axis). The y- 
axis is distance parallel to the line joining the Sun to the Galactic 
Centre. The grid scale is 400 pc for each small square. The range 
of distances is from 4000 to 12000 pc along the line of sight, and 
from —2000 to +2000 pc in the x-axis. The origin is at the Sun. 




Figure 11. The same plot as Figure ^ but for 5 = —22° 



al. (1999: Sect. 7, Fig. B) also show this asymmetry derived 
from TMGS data. 



6.3 The three-dimensional bulge 

The morphology of the bulge can be examined by fitting 
the isodensity surfaces to D{r} = D(r,l,h). The results of 



the previous subsection argued for non-axisymmetry in the 
bulge, so the next stage was to determine the parameters. 

Ellipsoids were used for the fit, with two axes in the 
Galactic plane and a third perpendicular to these. The pos- 
sible tilt of the bulge out of the plane was neglected as there 
is no evidence for this (Weiland et al. 1994). Also the posi- 
tion of the Sun 15 pc above the plane (Hammersley et al. 
1995) does not have a significant infiuence since the bulge 
extends much further from the plane. 

The Galactocentric distance along the major axis for 
different isodensity ellipsoids is 



(22) 



and the distance along the minor axis is t/Kz. The pro- 
jections of the vector distance to the Galactic centre are 
represented xi and X2 (Fig. ^2|), and z is the distance to the 
plane. K2 and are the axial ratios between axes xi and 
X2, and xi and 2, respectively. Both ratio are defined to be 
greater than one. 

From the same figure Xi and X2 are defined as follows: 



xi — Rcoa{P — a) 
and 

X2 — -Rsin(/3 — a), 
with 

R = (r cos b)^ + Rq — 2rRo cos b cos I, 
z = r sin b, 

and, following the sine rule, 



_i r cos 6 sin / 
R ' 



(23) 

(24) 

(25) 
(26) 

(27) 



The ellipsoids have four free parameters: Ro, the Sun- 
Galactic centre distance (the ellipsoids are then centred on 
this position); and Ky, the axial ratios with respect to 
the major axis (x); and a, the angle between the major axis 
of the triaxial bulge and the line of sight to the Galactic 
centre {a between 0° and 90° is where the tip of the major 
axis lies in the first quadrant). 

Three-dimensional ellipsoids are fitted to 20 isodensity 
surfaces (from 0.1 to 2.0 star pc~"^, in steps of 0.1) with the 
four free parameters. 

The four parameters are then averaged for the 20 ellip- 
soids and the results are: 

Ro = 7860 ± 90 pc, 

K2 = 1.87 ±0.18, 

Kz = 3.0 ±0.9 

and 



Q = 12 ± 6 deg. 



(28) 



The errors are calculated from the average of the ellip- 
soids and so do not include possible systematic errors (for 
example: subtraction of the disc, contamination from other 
components, methodological inaccuracies of the inversion, 
etc.), which are difficult to determine. However, by far the 
largest effect on the bulge counts is the massive asymme- 
try in the counts caused by the triaxiality of the bulge, as 
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Figure 12. Cut of an ellipsoidal bulge in the Galactic plane. C is the Galactic centre, P is a given point on the ellipsoid and S is the 
Sun. 



shown in Hammersley et al. (1999). The other systematic ef- 
fects are at least an order of magnitude below this, so while 
they do have an effect, it is small. Hence, the true errors are 
larger than stated but tests suggest that they do not alter 
the general findings presented here. 

The error in Kz is quite large and is due to the non- 
constant ajcial ratio of the ellipsoids. tends to increase 
towards the centre, i.e. the outer bulge is more circular than 
the inner bulge. This will be further discussed in §3.4. 



6.3.1 Axial ratios and orientation 

The axial ratios of the bulge are 1:0.54:0.33. These numbers 
indicate that the bulge is triaxial with the major axis close 
to the line of sight towards the Galactic centre. In general, 
the result presented here are in agreement with those from 
other authors. The projection, as viewed from the position 
of the Sun, of an ellipsoid of the above characteristics, gives 
an ellipse with axial ratio 1.7±0.5 (i.e. 1:0.58). This is com- 
patible with the value of 1:0.6 obtained by Weiland et al. 
(1994) or 1:0.61 by (Kent et al. 1991). 

From a dynamic model assuming a gas ring in a steady 
state, Vietri (1986) finds axial ratios of 1:0.7:0.4, which is 
close to our result. Binney et al. (1991) found q = 16 deg for 
a bar, i.e. a triaxial structure in the centre of the Galaxy, in 
order to explain the kinematics of the gas in the centre of the 
Galaxy. Weinberg (1992) gives a = 36 ± 10 deg and K2 = 
1.67 from his analysis of IRAS data.. More recently, Nikolaev 
& Weinberg (1997) obtained a bar from /iJylS' sources with 
a = 19 deg and K2 between 2.2 and 2.7. Stanek et al. (1997), 
based on the analysis of optical photometric data for regions 
of low extinction, predicted an a between 20 and 30 deg and 
1:0.43:0.29 axial ratios, which is also quite close to the result 
presented here. Various authors have examined the COBE- 
DIRBE flux maps for triaxiality: Dwek et al. (1995) give 
higher eccentricity values for the axial ratios, 1:0.33:0.22, 



but the angle a = 20 ± 10 deg is compatible with the value 
given here; Binney et al. (1997) derive 1:0.6:0.4 ratios and 
an angle a 20 deg; Freudenreich (1998) obtained a best fit 
with 1:0.38:0.26 ratios and a = 14 deg. Normally, when the 
low latitudes are excluded the fit of the triaxial bulge has 
an angle of about 25 degrees (Sevenster et al. 1999). The 
majority of the above authors did not use inversion; rather 
they fitted the fiux or the star counts to models using a 
priori assumptions. Binney et al. (1997) were the exception 
in using inversion on the C05i5-DIRBE surface-brightness 

map^ and, on the basis of specific assumption, obtained 
results close to those presented here. 



6.3.2 Galactocentric distance 

The distance Ro derived here is slightly less than that used 
in the W92 model of the disc (8.5 kpc). However, the small 
changes in Rq can be compensated by small changes in the 
other model parameters, such as the scale length, so that 
the predicted counts remain the same. As the model used 
already gave a good fit to the disc, we decided not to make 
ad hoc modifications to account for a smaller Ro since the 
disc is not the subject in this paper. 

The lack of previous assumptions makes the determi- 
nation of Ro presented here different from those of other 
authors. In particular, no information is required on the ob- 
jects observed. However, the values determined here are very 
close to the currently accepted value of just under 8 kpc. 
Reid et al. (1988) deduce a value Ro = 7.1 ± 1.5 kpc from 



II The inversion of the flux and the inversion of the star counts 
are significantly different. Since star counts provide a function for 
each region of space and the flux is only one number for each of 
those regions, the inversion of the flux is less suitable for directly 
extracting information from the data and further assumptions are 
needed. 
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Table 4. Relationship between the maximum distance of the el- 
lipsoid and the bulge star density. 



c 


n 




n 




IPC I 




VPC 1 


3020 ± 810 


0.1 


1620 ± 250 


1.1 


2630 ± 650 


0.2 


1580 ± 240 


1.2 


2420 ± 590 


0.3 


1540 ± 240 


1.3 


2230 ± 490 


0.4 


1460 ± 230 


1.4 


2120 ± 450 


0.5 


1420 ± 230 


1.5 


1990 ± 380 


0.6 


1390 ± 230 


1.6 


1900 ± 350 


0.7 


1380 ± 240 


1.7 


1840 ± 330 


0.8 


1360 ± 250 


1.8 


1720 ± 280 


0.9 


1360 ± 260 


1.9 


1670 ± 270 


1.0 


1320 ± 220 


2.0 



direct observations of Sgr B2. Gwinn et al. (1992), by means 
of observations of masers in W49, derive Rq = 8.1 ± 1.1 
kpc. Moran (1993) obtains Ro = 7.7 kpc, from OH/IR stars 
distances. Turbide & Moffat (1993) obtain Ro = 7.9 ± 1.0 
kpc, from measurements of the distances to young stars by 
means of CCD photometry and assuming that there is no 
metalhcity gradient in the outer regions of the Galaxy; al- 
though they get 7.2 kpc when a certain gradient is assumed. 
Paczyhski & Stanek (1998) derived Ro = 7.97 ± 0.08 (sys- 
tematic effects make the true error larger) from the compar- 
ison between Hipparcos and OGLE data. Oiling & Merrifield 
(1998a, 1998b) obtain Ro = 7.1 ±0.4. Etc. Generally, many 
studies based on indirect measurements claim the Galacto- 
centric distance to be somewhat less than 8.0 kpc (see also 
the review by Reid, 1993). 



6.3.3 Density as a function of the distance to the Galactic 
centre 

A power law with exponent —1.8 is observed in the cen- 
tre of the bulge and also in other galaxies (Becklin & 
Neugebauer 1968; Sanders & Lowinger 1972; Maihara et 
al. 1978; Bailey 1980; see review by Sellwood & Sanders 
1988). When the density function D{t) (Table ^ is fitted 
to D{t) = A{t/toy-^ exp{-{t/toy'), with 7, io and A as free 
parameters, then we obtain 

D{t) = 1.17(t/2180 pc)"^'^exp(-(t/2180 pc)^'^) 

starpc"^. (29) 

This gives an estimate of the fall-off in density between 
1.3 and 3.0 kpc from the centre in the direction parallel to 
the major axis or between 0.4 and 1.0 kpc in the direction 
perpendicular to the plane. As can be seen in Fig. |l^, the 
dispersion of points around this law is large, so it is possible 
to accommodate other functions or even a different set of 
parameter. A different luminosity function amplitude would 
change the amplitude of the stellar density, A. If the nor- 
malization for the luminosity function were incorrect then 
the factor needed to multiply the luminosity function would 
be used to divide the star-density amplitude. 



Ie+01 




0.0 1000.0 2000.0 3000.0 4000.0 

t(pc) 



Figure 13. Fit of the density distribution. The solid line is the 
best fit using eq. (H). 
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Figure 14. N{mK = 9.0 mag) along the three strips that are 
used with constant declinations: <5 = —30°, which cuts the plane 
at I = —1°; 5 = —22°, which cuts the plane at Z = 7°; and 
5 = —16°, which cuts the plane at Z = 15° once the W92 disc and 
bulge (according to eqs. (|28|) and (p9|)) are subtracted. 



6.3.4 Goodness of the inversion 

The residual counts for mA' < 9 after subtracting both the 
bulge determined here and the W92 disc model from the 
original counts are plotted in Figure |l^. As can be seen, 
the off-plane residual counts (the |6| < 2° regions are clearly 
contaminated by other components) are reduced to typically 
a few per cent of the original counts shown in Figure ^. For 
the 5 = —30° strip the residuals are typically 100 star/deg^ 
compared to the 1500 star/deg^ in the original counts. Hence 
the proposed bulge parameters do accurately reproduce the 
observed counts. 
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6.3.5 A triaxial bulge 

From Figs. ^ and the non-axisymmetry was determined 
for tiie plane. Furtliermore, tlie axial ratio K2 is close to 2 
(and not 1, the condition of axisymmetry) . Therefore the 
bulge is a triaxial ellipsoid orientated in such a way that the 
minor axis is perpendicular to the Galactic plane, and the 
angle between the major axis and the Sun-Galactic centre 
line is 12° in the first quadrant. 

Whether this structure is called a bar or triaxial bulge 
is not only a question of wording. Apart from the morphol- 
ogy, the population is also has to take into account: bulges 
are older than bars (Kuijken 1996), though both are older 



5.1) would 



than the disc. Precise calculations of the age (see \ 
be necessary to differentiate between them. However, there 
is evidence of another lengthened structure, a bar, (Ham- 
mersley et al. 1994; Calbet et al. 1996; Garzon et al. 1997) 
whose angle is ~ 75° in the first quadrant. This has major 
star formation regions at both extremes (towards I — 27° 
and I — —22°) and there is evidence for a preceding dust 
lane (Calbet et al. 1996). If this other component exists 
then the structure discussed in this paper must be called 
a "bulge", unless we are prepared to entertain the notion 
that the Galaxy has two bars. 



6.4 Bulge with variable ellipsoids 

A large error in is obtained when it is assumed constant, 
as was indicated in the previous subsection. Therefore, it is 
possible that the values are not constant, and so another 
dependence on the isodensity contours was tried. When the 
ellipsoids are fitted allowing a linearly variable K^, then 



= (1.66 ± 0.17) + (1.73 ± 0.14)L> 



(30) 



(where the units of D are star pc~^), whose weighted average 
is Kz = 3.0, as obtained in eq. (|2q). The other parameters 
{K2, Ro and a) remain nearly constant with respect to D. 

This variation of Kz is independent of the trial solution 
in the iteration process (see § [4.l[ ). A fourth iteration was 
performed for both the luminosity function and the density 
with the feed-back of the variable K^, and it could be seen 
that the same parameters are recovered again, within a l-cr 
error. Indeed, the xi-z ratio is 



Kz = (1.76 ± 0.32) + (1.70 ± 0.27)L>. 



(31) 



This linear dependence is valid in the density interval 
from 0.1 to 2.0 star pc~'^. For highest densities, is ex- 
pected to grow more slowly. can also be expressed as a 
function of t, although this dependence is non-linear. The 
fit to an exponential law is: 



Kz = (8.4± 1.7)exp 



(2000 ± 920) pc 



(32) 



and is valid for the range of distances, t, used here (see 
Fig. |l^) . Figure ^ shows the variation of eccentricity as a 
function of the density. 

With Kz so defined, the density, D, is given by Table ^ 
or Figure |l6[P|. 



** The density is expressed as a function of t/Kz, the distance 
along the z-axis, because the variation of Kz with t fluctuates too 
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Figure 15. Cut of the bulge in the xi—z plane when Kz is given 
by eq. (^^. The eUipses represent isodensity lines between 0.1 
and 2.0 star pc^'^, with intervals of 0.1. 



Table 5. Relationship between the distance along the minor axis, 
t/Kz and the density, when Kz is given by (p^. 



t/Kz 


D 


t/Kz 


D 


(pc) 


(star pc~^) 


(pc) 


(star pc"'') 


1400 ± 380 


0.1 


470 ± 60 


1.1 


1170 ± 290 


0.2 


440 ± 50 


1.2 


1010 ± 240 


0.3 


430 ± 40 


1.3 


880 ± 200 


0.4 


410 ± 40 


1.4 


780 ± 160 


0.5 


390 ± 40 


1.5 


700 ± 140 


0.6 


380 ± 30 


1.6 


640 ± 120 


0.7 


370 ± 30 


1.7 


570 ± 100 


0.8 


360 ± 30 


1.8 


540 ± 90 


0.9 


350 ± 30 


1.9 


510 ± 80 


1.0 


330 ± 20 


2.0 


The best 


fit to a 


law of type 


D{t/Kz) 



A{t/{Kzto))-^-^eM-{i:/{Kz xto))^) is: 
f t/Kz 



D(t/Kz 



star pc 



0.106 



\1820 pc 



exp 



t/Kz 
1820 pc 



(33) 



6.4-1 Goodness of the inversion 

The residual counts after both the bulge determined here 
with variable Kz and the W92 disc model are subtracted 
from the original counts for m^ < 9 are plotted in Figure |l^ 
As can be seen the residuals are now somewhat lower than 
when Kz is constant (Fig. [l4| ). Typically the residuals are 
now 50 to 100 star/deg^ and the maximum has fallen from 
300 star/deg^ with constant Kz to 200 star/deg^. Therefore, 



much. The ellipsoid size decreases when the density, D, increases; 
however, Kz increases with D, so the axis xi increases. Hence, 
the variation of D as a function of t is too sensitive to noise. 
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Figure 17. N{mK = 9-0 mag) along the three strips that are 
used with constant dechnations: & = —30°, which cuts the plane 
at / = —1°; <5 = —22°, which cuts the plane at Z = 7°; and 
5 = —16°, which cuts the plane at Z = 15° once the W92 disc and 
bulge (according to eqs. (|28|), and (^)) arc subtracted. 



the variable does provide a better fit to the observed 
counts. 

The aspect of the bulge as seen by an observer far away 
in the z-axis, i.e. the Galaxy observed face-on, would be 
as shown in Figure ^| The sharp fall-off in density is very 
noticeable. The bulge in a face-on Milky Way-like Galaxy 
presents, according to our results, a very high contrast be- 
tween central regions (with up to 10000 star pc^^) and re- 
gions at 3 kpc in the major axis (with 100 star pc~^). 

Whether this variation of Kz is a true feature of the 
density distribution or not is a matter for further investiga- 
tion. However, this is observed in other galaxies (Varela et 
al. 1993) and we do not believe that the result of this sub- 
section is due to systematic errors, although this possibility 
cannot be totally excluded. 

Other possible causes for this might be either that a 



Figure 18. Projection of the bulge, eq. (^), when it is observed 
face-on (integration of z direction). The square is 8 kpcx8 kpc 
centred on the Galactic Centre. The outer contour represents 10 
star pc~^, the second contour stands for 410 star pc~^, etc., and 



the inner contour stands for 4810 star pc 
consecutive contours is 400 star pc~'^). 



(the interval between 



superposition of two components is being observed, i.e the 
bulge and another structure, a bar, closer to the plane. If this 
were true, the luminosity function would have two different 
populations, especially in the regions closest to the plane. A 
gradient within the bulge is also possible. A greater number 
of bright stars in the innermost bulge (as observed by Calbet 
et al. 1995), with a smooth variation from the inner to the 
outer bulge, could be responsible of this effect. Both of these 
causes would lead to a gradient in the luminosity function. 
However, as the luminosity function has been assumed to 
be constant the result after inversion would be a gradient in 
Kz- Tests on the data indicate that this is possible. Giving 
the luminosity function a gradient in z, but such that the 
luminosity function remains within the error bars for a de- 
termined average function, is sufficient to produce changes 
in the observed gradient in Kz- In any case, the errors in 
the luminosity function (see Table ^ do limit this variation 
of populations. 



6.5 Stellar content of the bulge 

Integrating the density over all space will give the stellar 
content of the whole bulge: 



iV 



D{t)dV- 



(34) 



The volume element dV , under a change of variable to 
elliptic coordinates t, 9 and 0, is related to the Cartesian 
coordinates x = t sin 6 cos 4>, y ~ sin 9 sin cf), z = cos T 
through 

dV = hthgh,i,dtd9d(j>, 



(35) 
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+ 



dy 
dqi 



+ 



dqi 



(36) 



Hence, 



dV — t sm9\ / sin- 



+ cos^ 9 



cos^ (j) + 



/ cos^ 9 .On 



COS^ (j) + ■ 



dtd9d(j). 



(37) 



The result is 2.8 x 10^" stars for ^ 3 with D from 
eg. (p9|); and 4.1 x 10^" stars with a variable from eq. 
( P^ and D from eq. (js^), i.e. a factor 1.4 greater. This is, of 
course, only an estimation which includes an extrapolation 
of D to all space and the assumption of a correct luminosity- 
function normalization (see §P). Nevertheless, it leads to an 
order of magnitude for the mass of the Galactic bulge (taking 
an average mass for a star of ~ IMq) compatible with other 
data (for instance, ~ 2 x lO^'^ M© in Gould 1997); so this 
supports the normalization and the extrapolation. From the 
integration of the luminosity function it is found that TMGS 
stars from the whole bulge {niK < 9.0 mag) represent only a 
fraction (~ 2 x 10^^) of the total number of stars, i.e. 6 x 10^ 
stars for Kz ~ 3. 



7 HOW DIFFERENT WOULD THE RESULTS 
FOR A DIFFERENT DISC MODEL BE? 

Errors in different parts of the inversion procedure used here 
will lead to changes in the results. One important source of 
error may be the disc model that is used: were a different 
model to be used, the answers would be different. Clearly 
the answer to a certain extent depends on the new model 
As has been detailed earlier in 



3.3 



the disc model 



to use. 

used here is in good agreement with the observed TMGS 
star counts where the disc is isolated and so the expectation 
is that its extrapolation to regions where bulge and disc are 
observed will lead only to small errors in the counts. By 
definition, the bulge is an excess over the extrapolated disc 
in central regions of the Galaxy so, also by definition, the 
error of the present disc model cannot be very large once its 
fitting to observational data in external parts of the Galaxy 
has been tested. 

From the integral equation it can be deduced that 
these errors, SNK,disc{'mK), follow for all regions {l,b): 

SNK.disciniK) = -to 

n 2 

X / 5$if,buigc(mA- + 5 — 51ogi„Pif)Abuigc,K(pK)pj<:rfpK 



(38) 



-UJ I $K,bulge(mx + 5 - 51ogio pK)(5Abulgo,A'(pK) 

Jo 

X pjfdpK, 



If we know 5NK,disc{'mK), which differentiates the "real 
disc" from our model, we could derive how large (5$ and 
SA are. The inversion procedure explained in 94] produces 
solutions which are close when we begin the iteration from 
counts that are similar, as can be seen from eq. (|l^). Hence, 
for small SNK,diBc{'mK), 5$ and 5A are also small. That is, 
the behaviour is not a chaotic such that small departures 
from the original counts would not produce very different 
solutions. 

For instance, let us suppose that there is an error, SHr, 
in the scale length of the disc (equal to 3.5 kpc in the W92 
model we assumed). This leads to an error in the density 
due to the disc of 



(39) 



where R is distance from the Galactic centre and Rq is this 
distance for the Sun (8.5 kpc in the W92 model). Hence, by 
means of eqs. (^, (^ and (^) for the disc. 



SNK,diBc{mK) ^ to 



5hR 

h 



R Jo 



"I'K,disc(mK + 5 - 51og^(, pk) 



X Adiscif (p/f )(-R - RQ)pKdpK, 



(40) 



which can be set equal to expression ( |38| ) for all ttik, I and 
b. However, whilst in principle the change in the disc is pro- 
portional to the scale length, and there are certainly values 
quoted in the literature as low as 2.2 kpc (Ruphy et al. 
1996), it should be remembered that it is already known 
that the W92 model gives an excellent fit in the areas where 
the disc is isolated. Therefore, were one to alter the scale 
length, then other parameters also would have to be varied 
to compensate, otherwise the excellent agreement would be 
lost. It would be difficult to change the disc more than a few 
per cent without the effect becoming noticeable. 

In the selected regions the bulge counts are dominant. 
For instance, the maximum contribution of the disc in the 
region {I = 0.3°, b = -2.0°) is 1200 star/deg^ up to 
9th K-magnitude whereas the total counts are around 6000 
star/deg^ (Fig. i.e. in this case only 20% of the sources 
are from the disc. The ratio varies according to the region 
examined, but in most of the regions used the bulge is the 
dominant feature. Furthermore the error in the number of 
bulge sources is determined by the error in the number of 
disc sources, therefore if the relative proportion of the disc 
sources is low and the disc model gives a good fit to the 
TMGS counts this implies that the error introduced to the 
bulge counts will be of the order of a few per cent, probably 
below the Poissonian noise. Therefore, the errors in the disc 
affects the shape and luminosity function of the bulge only 
slightly. 



7.1 Experiments of inversion varying the 

parameters of the disc or the extinction 

A simple test can be carried out to verify what has been 
claimed in this section: small changes in the parameters of 
the disc (or also the extinction) do not greatly affect the 
results, i.e. there is no chaotic behaviour. 

We run the same inversion programs again to obtain 
both the luminosity function and the density distribution. 
Two examples are shown in this subsection: a) inversion with 
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Figure 19. Luminosity function with different parameters for the 
disc model and the extinction as well as for the W92 model of the 
bulge. Reference luminosity function is the one obtained in 

ha = 3.0 kpc instead of the original value of hn = 3.5 
kpc; b) inversion with extinction normalization coefficient 
Ak ~ 0.05 mag kpc~^ instead of Ak = 0.07 mag kpc"\ The 
new luminosity functions are shown in Fig. |l^ in comparison 
with that obtained in section |^. 

Both inversions with new disc and extinction are fit 
to constant axial-ratio triaxial ellipsoids respectively with 
parameters: a) Ro = 8400 ± 190 pc, = 2.5 ± 1.3, Ky = 
1.75 ± 0.05, a = 12° ± 3°; b) Ro = 7600 ± 130 pc, = 
4.1 ± 1.1, Ky = 1.70 ± 0.05^0 = 9° ± 2°. These values are 
close to those obtained in §pj, which is an indication of the 
robustness of the method of inversion. 

In case a), the luminosity function for very bright stars 
in K is higher than the reference one in comparison with 
the faintest parts due, perhaps, to a defect of outer bulge 
stars. The disc model in a) is unrealistic and provides fur- 
ther star counts in the Galactic centre than there should be 
(~ 25% more stars than in the reference model); the outer 
regions of the bulge would have zero, or negative, counts 
once the disc is subtracted, so they do not contribute to the 
weighted average of the luminosity function. In case b), the 
Galactic centre is closer to us, as expected if the extinction 
is lower. No physical meanings can be derived from these ex- 
periments, since the disc model in a) or the extinction model 
in b) is less exact than that in the reference case. They sim- 
ply provide a verification of the robustness of the inversion 
method. 



8 CONCLUSIONS 

The procedure used here is rather different from that of 
those authors who fit the parameters directly to the star 
counts. First, the counts were inverted. Then, after the lumi- 
nosity function and density distribution were evaluated and 
an approximate ellipsoidal shape was evident, the parame- 
ters could be fitted for each isodensity surface. Assuming an 
ellipsoidal bulge with constant parameters for all isodensity 
regions and fitting these parameters to the counts is less rig- 



orous since there is no a priori evidence for this assumption. 
In fact, our method suggests that constant parameters for 
the ellipsoids do not give the best fit for the density, D{r). 
Instead, a decreasing major-minor axial ratio from inside to 
outside would provide best results. 
These results are: 

The distance to the centre of the bulge, i.e. the centre of 
the Galaxy, is 7.86 ± 0.09 kpc (systematic effects make the 
true error larger). 

The relative abundance of the brightest sources in the 
bulge {Mk < — 8.0 mag) is much less than in the disc. 

The bulge is triaxial with axial ratios 1:0.54:0.33, the mi- 
nor axis perpendicular to the Galactic plane, and the major 
axis nearly along the line of sight to the Galactic centre. The 
best fit giving an angle equal to 12 deg shifted to positive 
Galactic longitudes in the plane in the first quadrant. 

A gradient in the major-minor axial ratio is measured. 
However, there are various possible caused which include 
eccentricity of the true density-ellipsoid gradient or a popu- 
lation gradient. 

The stellar density drops quickly with distance from 
the Galactic centre (i.e. the density distribution is sharply 
peaked). The —1.8 power-law observed at the Galactic cen- 
tre needs to be multiplied by an exponential to account for 
the fast drop in density in the outer bulge. 
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